Exercise 1: Finding Michel Electrons¶

Michel electrons come out of a muon decay. They are electrons, so they look like an electromagnetic shower coming out of a long muon track.

Note

Need refresh to make sure we get a decent spectrum at the end of the tutorial.

Requirements:¶

U-ResNet for Semantic Segmentation

I. Motivation¶

Michel electrons are used in LArTPC experiments as “candles”: we understand their energy spectrum very well, so they can be used to calibrate the detector and to compare the quality of different detectors. Michel electrons spectrum is one of the first “tokens” that an experiment will show to prove that the detector is running and we can successfully reconstruct the data. Since their energy is in a range up to ~50MeV, they are also representative of the detector response to electromagnetic particles generated by neutrino activity.

https://encrypted-tbn0.gstatic.com/images?q=tbn:ANd9GcQsezzzPeMfyra-aaqqooRh-IAVN61XbgDDuQ&usqp=CAU

Fig. 9 Example of Michel electron spectrum¶

What are we looking for?¶

Fig. 10 Energy loss per unit distance (MeV/cm) for electrons traveling in liquid argon. From [AAA+17].¶

The primary Michel electron that comes out of the muon decay-at-rest can lose energy either through ionization (collision stopping power), or through the production of radiative photons via Bremsstrahlung (radiative stopping power). Beyond a certain energy, the radiative stopping power is greater than the collision stopping power. If the radiative photons have enough energy, they can pair-produce, i.e. turn into an electron-positron pair. They, in turn, can produce new radiative photons, and so on. They can also undergo Compton scattering. A cascade of electrons and photons (electromagnetic shower) happens.

Ionization produces track-like energy depositions, whereas the photons can travel some distance before converting into a secondary electron. Hence Michel electrons have two clear topological features: a primary ionization, which is track-like at the end of the muon track, and some scattered energy deposits much further away which come from these radiative photons.

In this exercise, we are only concerned with finding the primary ionization of the Michel electron. We will purposefully ignore for now the radiative photons.

Fig. 11 Example of Michel electron topology¶

II. Setup¶

We first need to set the working environment and the path to the validation dataset.

import os, sys
SOFTWARE_DIR = '%s/lartpc_mlreco3d' % os.environ.get('HOME') # Path to your `lartpc_mlreco3d`
DATA_DIR = os.environ.get('DATA_DIR')
# Set software directory
sys.path.append(SOFTWARE_DIR)

We will also need Plotly for visualization:

import numpy as np
import matplotlib.pyplot as plt 
import yaml

import plotly
import plotly.graph_objs as go
from plotly.subplots import make_subplots
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
init_notebook_mode(connected=False)

The imports below load some auxiliary functions and classes required to run the full chain in interactive mode.

from mlreco.main_funcs import process_config, prepare
import warnings
warnings.filterwarnings('ignore')

/usr/local/lib/python3.8/dist-packages/MinkowskiEngine/__init__.py:36: UserWarning:

The environment variable `OMP_NUM_THREADS` not set. MinkowskiEngine will automatically set `OMP_NUM_THREADS=16`. If you want to set `OMP_NUM_THREADS` manually, please export it on the command line before running a python script. e.g. `export OMP_NUM_THREADS=12; python your_program.py`. It is recommended to set it below 24.

Let’s load the configuration file to setup the full chain architecture and weights.

cfg = yaml.load(open('%s/inference.cfg' % DATA_DIR, 'r').read().replace('DATA_DIR', DATA_DIR),Loader=yaml.Loader)
process_config(cfg, verbose=False)

Config processed at: Linux tur015 3.10.0-1160.42.2.el7.x86_64 #1 SMP Tue Sep 7 14:49:57 UTC 2021 x86_64 x86_64 x86_64 GNU/Linux

$CUDA_VISIBLE_DEVICES="0"

This cell loads the dataset and the model to the notebook environment. Optionally, you can specifiy the list of images by listing the image ID numbers:

prepare(cfg, event_list=[0, 1, 2, ...]

If the model’s trained weights have been loaded correctly, it should display a message like this ($PATH_TO_WEIGHTS will be replaced by wherever you placed your trained full chain weights):

Restoring weights for from $PATH_TO_WEIGHTS… Done.

# prepare function configures necessary "handlers"
hs = prepare(cfg)
dataset = hs.data_io_iter

Welcome to JupyROOT 6.22/09

Loading file: /sdf/home/l/ldomine/lartpc_mlreco3d_tutorials/book/data/mpvmpr_012022_test_small.root
Loading tree sparse3d_reco
Loading tree sparse3d_reco_chi2
Loading tree sparse3d_pcluster_semantics_ghost
Loading tree cluster3d_pcluster
Loading tree particle_pcluster
Loading tree particle_mpv
Loading tree sparse3d_pcluster_semantics
Loading tree sparse3d_pcluster
Loading tree particle_corrected

Warning in <TClass::Init>: no dictionary for class larcv::EventNeutrino is available
Warning in <TClass::Init>: no dictionary for class larcv::NeutrinoSet is available
Warning in <TClass::Init>: no dictionary for class larcv::Neutrino is available

        Since one of the GNNs are turned on, process_fragments is turned ON.
        
        Fragment processing is turned ON. When training CNN models from
         scratch, we recommend turning fragment processing OFF as without
         reliable segmentation and/or cnn clustering outputs this could take
         prohibitively large training iterations.
        
        Since one of the GNNs are turned on, process_fragments is turned ON.
        
        Fragment processing is turned ON. When training CNN models from
         scratch, we recommend turning fragment processing OFF as without
         reliable segmentation and/or cnn clustering outputs this could take
         prohibitively large training iterations.

Ghost Masking is enabled for UResNet Segmentation
Ghost Masking is enabled for MinkPPN.

Restoring weights for  from /sdf/home/l/ldomine/lartpc_mlreco3d_tutorials/book/data/weights_full_mpvmpr_012022.ckpt...

Done.

We can perform one forward iteration of the chain as follows:

data, result = hs.trainer.forward(dataset)

Segmentation Accuracy: 0.9824
PPN Accuracy: 0.8162
Clustering Accuracy: 0.9279
Shower fragment clustering accuracy: 0.9255
Shower primary prediction accuracy: 0.6000
Track fragment clustering accuracy: 0.9480
Interaction grouping accuracy: 0.9709
Particle ID accuracy: 0.8167
Primary particle score accuracy: 0.8864

The input data and label information are loaded onto the data_blob variable, while the outputs from the chain are stored inside result.

III. Setup the Evaluator¶

For the tutorial, we will use the analysis tools developed for easier handling of the full chain outputs.

from analysis.classes.ui import FullChainEvaluator

As we did before in the Analysis Tools tutorial, let’s setup the FullChainEvaluator with deghosting turned on:

predictor = FullChainEvaluator(data, result, cfg, deghosting=True)

entry = 8

IV. Selecting Michel Electrons¶

Step 1: Extract the semantic segmentation predictions.¶

Let’s first import plotting libraries for visualization:

from mlreco.visualization import scatter_points, scatter_cubes, plotly_layout3d
from mlreco.visualization.plotly_layouts import white_layout, trace_particles, trace_interactions, dualplot

The evaluator will handle ghost masking internally, so all you have to do is retrieve the predicted and true semantic segmentation labels for visualization and analysis:

pred = predictor._fit_predict_semantics(entry)
truth = predictor.get_true_label(entry, 'segment_label')
points = predictor.data_blob['input_data'][entry]
# Check if dimensions agree
assert (pred.shape[0] == truth.shape[0])

Let’s plot the true and predicted semantic labels side-by-side:

trace1, trace2 = [], []
edep = points[:, 5]

trace1 += scatter_points(points,
                        markersize=1,
                        color=truth, 
                        cmin=0, cmax=5, colorscale='rainbow')

trace2 += scatter_points(points,
                        markersize=1,
                        color=pred, 
                        cmin=0, cmax=5, colorscale='rainbow')

fig = dualplot(trace1, trace2, titles=['True semantic labels (true no-ghost mask)', 'Predicted semantic labels (predicted no-ghost mask)' ])

iplot(fig)

The true labels are plotted in the left, predicted on the right

Step 2. Identify muons and electrons.¶

By default, the label for tracks and michel electrons are 1 and 2, respectively.

track_label = 1
michel_label = 2

particles = predictor.get_particles(entry, primaries=False)
true_particles = predictor.get_true_particles(entry, primaries=False, verbose=False)

Let’s look at the list of reconstructed (predicted) particles:

trace1 = trace_particles(particles)
trace2 = trace_particles(true_particles)

fig = make_subplots(rows=1, cols=2,
                    specs=[[{'type': 'scatter3d'}, {'type': 'scatter3d'}]],
                    horizontal_spacing=0.05, vertical_spacing=0.04)
fig.add_traces(trace1, rows=[1] * len(trace1), cols=[1] * len(trace1))
fig.add_traces(trace2, rows=[1] * len(trace2), cols=[2] * len(trace2))
fig.layout = white_layout()
fig.update_layout(showlegend=False,
                  legend=dict(xanchor="left"),
                 autosize=True,
                 height=600,
                 width=1500,
                 margin=dict(r=20, l=20, b=20, t=20))
iplot(fig)

And the true particles:

true_particles

[TruthParticle( Batch=8   | ID=0   | Semantic_type: Track           | PID: Muon     | Primary: 1  | Interaction ID: 3  | Size: 697   ),
 TruthParticle( Batch=8   | ID=1   | Semantic_type: Track           | PID: Proton   | Primary: 1  | Interaction ID: 3  | Size: 222   ),
 TruthParticle( Batch=8   | ID=2   | Semantic_type: Shower Fragment | PID: Photon   | Primary: 1  | Interaction ID: 3  | Size: 904   ),
 TruthParticle( Batch=8   | ID=3   | Semantic_type: Shower Fragment | PID: Photon   | Primary: 1  | Interaction ID: 3  | Size: 749   ),
 TruthParticle( Batch=8   | ID=19  | Semantic_type: Track           | PID: Muon     | Primary: 1  | Interaction ID: 2  | Size: 100   ),
 TruthParticle( Batch=8   | ID=20  | Semantic_type: Track           | PID: Muon     | Primary: 1  | Interaction ID: 4  | Size: 321   ),
 TruthParticle( Batch=8   | ID=21  | Semantic_type: Track           | PID: Muon     | Primary: 1  | Interaction ID: 1  | Size: 2159  ),
 TruthParticle( Batch=8   | ID=24  | Semantic_type: Delta Ray       | PID: Electron | Primary: 0  | Interaction ID: 1  | Size: 23    ),
 TruthParticle( Batch=8   | ID=25  | Semantic_type: Michel Electron | PID: Electron | Primary: 0  | Interaction ID: 4  | Size: 85    )]

Ok, so we see that there is one predicted and true Michel electron in this image. We should first check if the two michel electrons actually correspond to one another. To see this, we first have to match all particles with one another:

matches = predictor.match_particles(entry, primaries=False)
def get_michels(matches):
    michels = []
    # Check if the predicted and true michels actually match to one another:
    for m in matches:
        if m[0].semantic_type == michel_label:
            michels.append(m[0])
    return michels

michels = get_michels(matches)

Step 3. Check if the Michel electron is adjacent to the end of a muon track.¶

To check if our candidate michel is actually adjacent to a track, we can set a minimum distance threshold between the michel electron and the track.

Wait, what if it was a misclassified Delta ray electron?

Delta ray electrons are knock off electrons that can happen along the trajectory of a muon, so if UResNet mispredicted the delta ray voxels as Michel voxels we would be wrong! To avoid that, let’s also make sure that the point of contact is at the end of the track. Again, many ways to do this check, this is just one possible heuristic.

If a Particle has a track semantic type, it will contain a (2 x 5) endpoints array as its attribute. The first three columns of endpoints are the voxel coordinates of the track’s endpoints. Optionally, we can use this information to be more selective.

I’ve wrote a function that checks if a michel electron is adjacent to a track in check_michel_adjacency as an example.

from scipy.spatial.distance import cdist

def check_michel_adjacency(michel_particles, particles, distance_threshold=10, check_endpoints=False):
    '''
    Given list of Michel electrons and list of particles,
    return a list of particle ids that are adjacent. 
    
    Inputs:
        - michel_particles (List[Particle]): list of michel electrons
        - particles (List[Particle]): list of particles 
        - distance_threshold: maximum distance to be considered adjacent
        - check_endpoints: option to use track endpoint information
        rather than Hausdorff distance.
        
    Returns:
        - adjacency_checks: list of particle ids that are
        adjacent to michels in <michel_particles> in order.
        If no adjacent tracks could be found, it will return 
        None instead of a particle id number (int). 
    '''
    # Check input michels are actually michels
    adjacency_checks = []
    for mp in michel_particles:
        is_adjacent = False
        if mp.semantic_type != 2:
            raise ValueError("Found non-michel particle {} with "\
                             "semantic type {} in list of michel electrons.".format(mp.id, mp.semantic_type))
        min_dist, adj_pid = np.inf, -1
        for p in particles:
            if p.semantic_type == 1:
                if check_endpoints:
                    dist = cdist(mp.points, p.endpoints[:, :3])
                else:
                    dist = cdist(mp.points, p.points)
                if dist.min() < min_dist:
                    min_dist = dist.min()
                    adj_pid = p.id
        if min_dist < distance_threshold:
            adjacency_checks.append(adj_pid)
        else:
            print("All candidate tracks are more than {}px apart from "\
                  "michel {}, minimum of {:.2f}px achieved for track {}".format(
                distance_threshold, mp.id, min_dist, adj_pid))
            adjacency_checks.append(None)
    return adjacency_checks

check_michel_adjacency(michels, particles, distance_threshold=10)

[4]

Step 4: Make a plot¶

Let’s use the voxel count as a substitute for the reconstructed energy. This approximation is fairly accurate for shower-like particles. Make a histogram with Michel electron candidates voxel counts.

Since there are only so many Michel electrons per entry, you will need to loop over several entries, possibly more than a batch size worth of entries. Let’s collect all michel electrons that have an adjacent track for 10 iterations-worth of data.

iterations = 10

all_candidates = []
for iteration in range(iterations):
    print("Iteration: {}...".format(iteration))
    data, result = hs.trainer.forward(dataset)
    # Initialize Evaluator"
    predictor = FullChainEvaluator(data, result, cfg, deghosting=True)
    for entry, index in enumerate(predictor.index):
        matches = predictor.match_particles(entry, primaries=False)
        particles = predictor.get_particles(entry, primaries=False)
        michels = get_michels(matches)
        adjacency = check_michel_adjacency(michels, particles)
        for i, m in enumerate(michels):
            if adjacency[i] is not None:
                all_candidates.append(m)

Iteration: 0...

Segmentation Accuracy: 0.9851
PPN Accuracy: 0.8156
Clustering Accuracy: 0.9136
Shower fragment clustering accuracy: 0.9040
Shower primary prediction accuracy: 1.0000
Track fragment clustering accuracy: 0.9334
Interaction grouping accuracy: 0.9463
Particle ID accuracy: 0.8030
Primary particle score accuracy: 0.9153

Particle 1 has no PPN candidates!
Particle 1 has no PPN candidates!

All candidate tracks are more than 10px apart from michel 24, minimum of 12.77px achieved for track 13
Iteration: 1...

Segmentation Accuracy: 0.9688
PPN Accuracy: 0.7979
Clustering Accuracy: 0.8609
Shower fragment clustering accuracy: 0.9102
Shower primary prediction accuracy: 0.7778
Track fragment clustering accuracy: 0.9569
Interaction grouping accuracy: 0.8940
Particle ID accuracy: 0.7564
Primary particle score accuracy: 0.8759

All candidate tracks are more than 10px apart from michel 19, minimum of 25.02px achieved for track 10

Iteration: 2...

Segmentation Accuracy: 0.9932
PPN Accuracy: 0.8321
Clustering Accuracy: 0.9033
Shower fragment clustering accuracy: 0.9487
Shower primary prediction accuracy: 1.0000
Track fragment clustering accuracy: 0.9720
Interaction grouping accuracy: 0.9106
Particle ID accuracy: 0.8269
Primary particle score accuracy: 0.9273

All candidate tracks are more than 10px apart from michel 13, minimum of 13.19px achieved for track 5
Particle 5 has no PPN candidates!
Particle 5 has no PPN candidates!
Iteration: 3...

Segmentation Accuracy: 0.9733
PPN Accuracy: 0.7615
Clustering Accuracy: 0.8377
Shower fragment clustering accuracy: 0.9434
Shower primary prediction accuracy: 0.9000
Track fragment clustering accuracy: 0.9534
Interaction grouping accuracy: 0.9414
Particle ID accuracy: 0.7778
Primary particle score accuracy: 0.9462

Particle 1 has no PPN candidates!
Particle 1 has no PPN candidates!
Particle 3 has no PPN candidates!
Particle 3 has no PPN candidates!
All candidate tracks are more than 10px apart from michel 16, minimum of 15.00px achieved for track 14
Particle 3 has no PPN candidates!
Particle 3 has no PPN candidates!
All candidate tracks are more than 10px apart from michel 9, minimum of 13.75px achieved for track 6

Iteration: 4...

Segmentation Accuracy: 0.9814
PPN Accuracy: 0.7792
Clustering Accuracy: 0.9067
Shower fragment clustering accuracy: 0.9516
Shower primary prediction accuracy: 0.5455
Track fragment clustering accuracy: 0.9500
Interaction grouping accuracy: 0.9290
Particle ID accuracy: 0.8947
Primary particle score accuracy: 0.9154

Particle 0 has no PPN candidates!
Particle 1 has no PPN candidates!
Particle 0 has no PPN candidates!
Particle 1 has no PPN candidates!

All candidate tracks are more than 10px apart from michel 12, minimum of 13.04px achieved for track 9
Iteration: 5...

Segmentation Accuracy: 0.9819
PPN Accuracy: 0.8271
Clustering Accuracy: 0.8859
Shower fragment clustering accuracy: 0.8787
Shower primary prediction accuracy: 1.0000
Track fragment clustering accuracy: 0.9530
Interaction grouping accuracy: 0.9193
Particle ID accuracy: 0.8868
Primary particle score accuracy: 0.9256

Particle 3 has no PPN candidates!
Particle 3 has no PPN candidates!
Particle 1 has no PPN candidates!
Particle 1 has no PPN candidates!
All candidate tracks are more than 10px apart from michel 16, minimum of 157.88px achieved for track 9

Particle 0 has no PPN candidates!
Particle 0 has no PPN candidates!
Iteration: 6...

Segmentation Accuracy: 0.9840
PPN Accuracy: 0.8182
Clustering Accuracy: 0.8362
Shower fragment clustering accuracy: 0.8671
Shower primary prediction accuracy: 0.9333
Track fragment clustering accuracy: 0.9604
Interaction grouping accuracy: 0.9478
Particle ID accuracy: 0.9091
Primary particle score accuracy: 0.8868

Particle 0 has no PPN candidates!
Particle 0 has no PPN candidates!

All candidate tracks are more than 10px apart from michel 11, minimum of 23.41px achieved for track 6

All candidate tracks are more than 10px apart from michel 8, minimum of 13.64px achieved for track 4
Iteration: 7...

Segmentation Accuracy: 0.9723
PPN Accuracy: 0.8146
Clustering Accuracy: 0.8628
Shower fragment clustering accuracy: 0.7941
Shower primary prediction accuracy: 1.0000
Track fragment clustering accuracy: 0.8976
Interaction grouping accuracy: 0.9419
Particle ID accuracy: 0.6957
Primary particle score accuracy: 0.8119

Particle 1 has no PPN candidates!
Particle 3 has no PPN candidates!
Particle 1 has no PPN candidates!
Particle 3 has no PPN candidates!
All candidate tracks are more than 10px apart from michel 6, minimum of 17.83px achieved for track 1
All candidate tracks are more than 10px apart from michel 7, minimum of 11.79px achieved for track 1
Particle 1 has no PPN candidates!

Particle 1 has no PPN candidates!
All candidate tracks are more than 10px apart from michel 9, minimum of 15.43px achieved for track 7
Iteration: 8...

Segmentation Accuracy: 0.9494
PPN Accuracy: 0.8093
Clustering Accuracy: 0.8473
Shower fragment clustering accuracy: 0.8639
Shower primary prediction accuracy: 0.6364
Track fragment clustering accuracy: 0.9450
Interaction grouping accuracy: 0.9534
Particle ID accuracy: 0.8000
Primary particle score accuracy: 0.9008

All candidate tracks are more than 10px apart from michel 14, minimum of 12.04px achieved for track 5
Particle 0 has no PPN candidates!
Particle 1 has no PPN candidates!
Particle 3 has no PPN candidates!

Particle 0 has no PPN candidates!
Particle 1 has no PPN candidates!
Particle 3 has no PPN candidates!
All candidate tracks are more than 10px apart from michel 15, minimum of 24.60px achieved for track 10
Iteration: 9...

Segmentation Accuracy: 1.0000
PPN Accuracy: 0.6214
Clustering Accuracy: 0.9372
Shower fragment clustering accuracy: 1.0000
Shower primary prediction accuracy: 0.0000
Track fragment clustering accuracy: 1.0000
Interaction grouping accuracy: 0.9000
Particle ID accuracy: 1.0000
Primary particle score accuracy: 0.8000

michel_voxel_counts = [p.points.shape[0] for p in all_candidates]

import matplotlib.pyplot as plt
import seaborn
seaborn.set(rc={
    'figure.figsize':(15, 10),
})
seaborn.set_context('talk')

plt.hist(michel_voxel_counts)
plt.xlabel("Voxel count")
plt.ylabel("Michel electron candidates")

Text(0, 0.5, 'Michel electron candidates')

Optional: how well did we do?¶

You can keep the exercise going by looking at the true Michel candidates (same heuristic, using the true semantic labels), matching them to the predicted Michel candidates and computing some metrics such as purity (fraction of predicted candidates that are matched) or efficiency (fraction of true Michel that find a match).

lartpc_mlreco3d Tutorials